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Abstract 

An  approximate  Riemann  solver  is  developed  for  the  governing  equations  of  ideal 
magnetohydrodynamics  (MHD).  The  Riemann  solver  has  an  eight-wave  structure, 
where  seven  of  the  waves  are  those  used  in  previous  work  on  upwind  schemes  for 
MHD,  and  the  eighth  wave  is  related  to  the  divergence  of  the  magnetic  field.  The 
structure  of  the  eighth  wave  is  not  immediately  obvious  from  the  governing  equations 
as  they  are  usually  written,  but  arises  from  a  modification  of  the  equations  that  is  pre¬ 
sented  in  this  paper.  The  addition  of  the  eighth  ave  allows  multi-dimensional  MHD 
problems  to  be  solved  without  the  use  of  sts^ered  grids  or  a  projection  scheme,  one  or 
the  other  of  which  was  necessary  in  previous  work  on  upwind  schemes  for  MHD.  A  test 
problem  made  up  of  a  shock  tube  with  rotated  initial  conditions  is  solved  to  show  that 
the  two-dimensional  code  yidds  answers  consistent  with  the  one-dimensional  methods 
developed  previously. 
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tract  No.  NASl-19480  while  the  author  was  in  residence  at  the  Institute  for  Computer  Applications  in  Science 
and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA  23681-0001. 
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1  Introduction 


The  governing  equations  of  ideal  magnetohydrodynamics  (MHD)  describe  the  physics  of  a 
conducting  fluid  in  which  the  following  assumptions  hold: 


where  p,  V,  r  and  L  are,  respectively,  characteristic  density,  speed,  time  and  length  scales 
for  the  problem,  c  is  the  speed  of  light,  and  c  and  a  represent  the  dielectric  constant  and 
conductivity  of  the  fluid.  These  equations,  written  in  conservation-law  form,  are 


puu  + 1  (p  -I-  -  BB 


^  B  uB  -  Bu 

\  j  [  (e  +  p+^)u-B{u-B)  ) 

where  I  is  a  3  x  3  identity  matrix,  p  is  the  density,  u  is  the  velocity,  p  is  the  pressure,  B  is 
the  magnetic  field,  and  E  is  the  energy,  defined  as 

p  u  •  u  B  •  B 

“  7-1  ^  2  2  ' 

Solutions  of  these  equations  can  yield  insight  into  a  number  of  problems  governed  by  fluid- 
dynamic  and  electromagnetic  effects. 

Much  of  the  past  work  in  solving  these  equations  has  b^n  based  on  Rusanov  and  Lax- 
Wendroff  techniques.  Only  recently  have  authors  begun  to  work  on  upwind  schemes  for 
solving  these  equations.  In  particular.  Brio  and  Wu  [2],  Zachary  and  Colella  [?],  and  Dai 
and  Woodward  [7]  have  done  some  of  the  early  development  of  Riemann-solver-based  schemes 
for  the  MHD  equations.  Their  work  has  been  based  not  on  the  system  of  eight  conservation 
laws  as  written  in  Equation  2,  but  instead  on  the  closely  related  system  that  comes  from 
assuming  =  constant  and  dropping  the  evolution  equation  for  Bx.  This  yields  a  7  x  7 
system.  The  reason  for  their  use  of  this  modified  system  arises  from  the  feM:t  that  one  of  the 
equations  governing  the  magnetic  field  is 

V  •  B  =  0  ,  (4) 
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Figure  1:  Waves  in  the  One- Dimensional  MHD  Riemann  Problem 


which,  in  one  dimension,  becomes  the  constridnt  Bg  —  constant. 

The  eigenvalues  and  eigenvectors  of  this  7x7  system  are  well  known  (see,  for  example, 
the  book  by  Jeffrey  and  Taniuti  [3]);  they  correspond  to: 

•  one  entropy  wave  traveling  with  speed  u; 

•  two  Alfven  waves  traveling  with  speed  u  ±  Co  where 


Bg 


is  the  Alfven  speed; 

•  four  magneto-acoustic  waves,  two  “fast”  and  two  “slow”,  traveling  with  speed  «  ±  c/ 
and  «  ±  c,  respectively,  where 


An  (x,t)  diagram  of  the  wave  interactions  at  a  cell  interface  is  shown  in  Figure  1. 

Given  these  seven  eigenvalues  and  corresponding  right  and  left  eigenvectors,  it  is  pos¬ 
sible  to  develop  a  linear  approximate  Riemann  solver  ala  Roe  [2,  ?],  or  a  more  nonlinear 
approximate  Riemann  solver  [7].  Once  some  questions  as  to  how  to  scale  the  left  and  right 
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eigenvectors  of  the  system  are  answered  (for  a  very  nice  description  of  the  problems  of  scaling 
in  the  MHD  eigensystem,  and  an  elegant  solution  to  these  problems,  see  the  paper  by  Rx)e 
and  Balsara  [5]),  a  robust  solver  for  one-dimensional  unsteady  problems  in  MHD  can  be 
developed. 

Building  a  code  capable  of  solving  two-  or  three-dimensional  problems  from  the  one- 
dimensional  Riemann  solver  building  block  is  not,  unfortunately,  as  straightforward  as  in 
the  case  of  the  Euler  equations.  In  the  one-dimensional  problem,  no  evolution  equation 
is  necessary  for  the  component  of  B  normal  to  a  cell  interface,  because  the  condition  of 
Equation  4  implies  that  8^  =  Bnn-  However,  in  a  two-dimensional  problem,  this  is  no 
longer  true.  In  two  dimensions,  the  discrete  constrsunt  corresponding  to  Elquation  4  is 

5;B«ds  =  0,  (5) 

facet 

and  so  a  jump  in  Bn  is  allowed  across  a  face;  it  simply  must  be  balanced  by  the  jumps 
across  the  other  faces  of  the  cell.  Thus,  a  separate  procedure  for  updating  this  portion  of 
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the  magnetic  field  must  be  implemented,  and  must  be  implemented  in  such  a  way  as  to 
satisfy  the  constraint  implied  in  Equation  5.  It  should  be  noted  that  the  V  •  B  constraint  is  a 
headache  not  just  for  upwind  schemes  for  MHD,  but  for  solution  of  MHD  problems  in  more 
than  one  dimension  by  any  method.  Typically,  one  of  three  approaches  is  taken  to  satisfy 
this  constraint: 

•  a  projection  scheme,  in  which  a  Poisson  equation  must  be  solved  to  subtract  off  the 
portion  of  the  magnetic  field  that  leads  to  a  non-zero  divergence; 

•  non-collocated  variables  (e.g.  a  staggered-grid  approach),  so  that  the  constraint  is  met 
identically; 

•  a  vector-potential  description  of  the  magnetic  field,  so  that  the  constraint  is  met  iden¬ 
tically. 

A  very  different  approach  is  taken  in  the  work  presented  here.  Instead  of  solving  a  seven- 
wave  Riemann  problem,  with  an  added  procedure  to  update  the  remaining  JB-field  component 
which  assures  that  Equation  4  is  satisfied,  ah  eight-wave  Riemann  solver,  in  which  all  of  the 
magnetic  field  components  are  updated,  is  developed  and  tested. 

2  Derivation  of  the  Eight-Wave  Riemann  Solver 

CJiven  the  primitive  variables 

W  =  (p,  u,  V,  w,  S*,  By,  Bz,p)  ,  (6) 
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Equation  2  may  be  rewritten  in  quasilinear  form  as 


dW 

-f-  Ap 

aw 

_  aw  _ 

aw 

0, 

(7) 
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The  Riemann  solver  would  normally  be  based  on  the  eigensystem  of  Ap, 

but  it  is  evident 

that  this  matrix  is  singular  — 

the  fifth  row 

of  the  matrix 

is  zero 

,  leading  to 

a  zero  eigenvalue. 

This  zero  eigenvalue  is  clearly  non-physical  —  the  eigenvalues  should  appear  either  singly  as 
the  z— component  of  the  flow  speed,  u,  or  in  pairs  symmetric  about  u.  It  also  does  not  bode 
well  numerically  —  the  mode  corresponding  to  this  eigenvalue  will  be  undamped. 

The  approach  taken  here  is  to  look  for  a  way  in  which  to  modify  the  governing  equations 
so  as  to  make  Ap  non-singular.  The  criteria  that  should  be  met  by  the  modified  matrix  >lp 
are: 

•  The  eigenvalues  and  eigenvectors  corresponding  to  the  seven  waves  in  the  one-dimensional 
(Bx  =  constant)  Riemann  solver  remain  unchanged; 

•  The  eigenvalue  corresponding  to  the  new  eighth  wave  is  equal  to  u  (the  only  physical 
choice  for  a  single  eigenvalue); 

•  The  left  and  right  eigenvectors  corresponding  to  the  new  eight  wave  “make  sense”; 

•  In  the  case  B*  =  constant,  the  eight-wave  Riemann  problem  reduces  to  the  seven-wave 
Riemann  problem. 

With  these  criteria  in  mind,  it  becomes  possible  to  find  a  modified  version  of  y4p,  given 
some  patience  and  some  facility  with  Maple’s  symbolic  manipulation  capabilities.  The  mod- 
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ified  matrix  that  meets  the  above  criteria 

is 
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The  eigensystem  of  this  matrix  is  composed  of  the  following  eight  waves,  with  their  corre¬ 
sponding  eigenvalues  A,  left  eigenvectors  i  and  right  eigenvectors  r: 


One  Entropy  Wave 


Two  Alfven  Waves 


A„  = 


L  = 


r»  = 


Ae  =  W 

K  =  (1,0,0, 0,0,0, 0,-1) 

fe  =:  (1,0,0,0,0,0,0,0)'^  ; 


u± 

y/P 

(0,  0,  — 5*,  By,0,  0)  ; 


Four  Magneto-acoustic  Waves 

A/,,  =  u±  c/,. 


—  (n  —  ^xBypCf,,  BxBzPCf^s  „  ^yP^f,a  ^2P^f,a 

f  -  („  B.B^Cjj  BiB.c/,,  Byp^j,  B,pc^,  \ 


(10) 


(11) 


(12) 
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One  “Divergence”  Wave 


Xj  =  u 

h  =  (0,0, 0,0, 1,0, 0,0) 

H  =  (0,0,0,0,1,0,0,0)'^  .  (13) 


It  is  important  to  note  that  the  first  seven  waves  yield  the  same  eigenvectors  and  eigen¬ 
values  as  the  seven-wave  Riem2um  problem,  with  the  additional  information  that  none  of 
them  carries  a  change  in  B*  (the  fifth  entry  of  each  right  eigenvector  is  zero),  and  none  of 
the  wave  strengths  is  proportional  to  a  jump  in  B,  (the  fifth  entry  of  each  left  eigenvector  is 
zero).  The  new  eighth  wave  travels  with  the  x— component  of  the  flow  speed  (its  eigenvalue 
is  u),  zmd  it  carries  a  jump  in  B»  (the  only  non-zero  entry  in  the  left  eigenvector  is  the  entry 
corresponding  to  B*),  and  affects  only  the  x— component  of  the  magnetic  field  (the  only 
non-zero  entry  in  the  right  eigenvector  is  the  entry  corresponding  to  Bx). 

It  is  clear  that  the  eigensystem  of  this  modified  matrix  has  all  of  the  desired  properties. 
In  the  case  B*  =  constant,  the  strength  of  the  eighth  wave  is  zero,  and  the  model  reverts  to 
that  of  the  seven-wave  problem.  The  new  wave  simply  gives  a  rational  procedure  for  dealing 
with  non-zero  jumps  in  B,  across  the  cell  faces,  which  will  in  general  occur  when  problems 
in  two  or  three  dimensions  are  being  solved.  The  question  remains,  however,  of  what  the 
modification  of  the  matrix  Ap  (and  the  corresponding  changes  to  Bp  and  Cp)  has  done  to 
the  system  of  conservation  laws. 

This  can  be  seen  by  collecting  the  source  terms  due  to  the  modifications  to  Bp  and 
Cp  and  transforming  to  conserved  variables.  The  new  equation  set,  which  heis  the  eight-wave 
eigensystem  described  above,  is 


(  \ 

'  0  ' 

P 

pu 

d 

p\x 

+  v. 

puu  -1- 1  (p  +  —  BnB 

B 

dt 

B 

uB  -Bu 

u 

<  ^  > 

^  (£;-l-p-l-^)u-B(u-B)  ; 

B^ 

(14) 


This  is  a  noteworthy  result:  the  source  term  that  must  be  added  to  Equation  2  is  proportional 
to  ^  'IB.  At  the  partial  differential  equation  level,  only  terms  that  are  equal  to  zero  have 
been  added  to  the  conservative  form  of  the  governing  equations.  So,  while  technically  the 
equations  are  no  longer  in  conservative  form,  the  deviations  from  conservation  will  be  very 
small.  It  is  only  by  writing  the  equations  in  this  slightly  non-conservative  form  that  the 
singularity  related  to  V  •  B  can  be  removed.  It  has  been  previously  noted  that  solving  the 
momentum  equation  in  non-conservative  form  can  remove  instabilities  related  to  non-zero 
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V  •  B  [1];  the  current  work  hopefully  reinforces  this  earlier  result,  and  sheds  further  light  on 
the  mechanism  for  stabilizing  the  equations,  as  well  as  applying  the  idea  in  a  novel  way  to 
develop  a  Riemann  solver  for  multi-dimensional  MHD. 

It  is  interesting  to  note  another  justification  of  this  particular  choice  of  source  term. 
Rewriting  Equation  2  slightly  by  expanding  some  of  the  terms,  the  following  form  of  the 
equations 


~  (pu)  -f  V  •  (puu)  +  v(p+  5^  j 


|  +  V.(pu)  =  0 


B  VB-BVjJB  =  0 


VB-t-BV  u-B- Vu-uV  B  =  0 
at 

B.V(uB)-(uB)VB  =  0  (15) 

is  obtained.  The  terms  that  are  proportional  to  V  •  B  have  been  underlined;  they  are  exactly 
the  same  as  the  source  term  defined  above.  Thus  it  can  be  seen  that  the  addition  of  the 
source  term  in  Equation  14  simply  acts  to  remove  the  terms  proportional  to  V  •  B  that 
appear  in  Equation  2. 

Another  interesting  note  is  what  the  evolution  equation  for  V  •  B  is  for  the  two  forms  of 
the  governing  equations.  This  may  be  seen  by  taking  the  divergence  of  the  evolution  equation 
for  the  magnetic  field  in  Equations  2  and  14.  For  the  original  form  of  the  equations,  the 
evolution  equation  is 


(dB 


VB  -b  BV  •  u  —  B  •  Vu  —  uV  •  B 


)  = 


^(V-B)  =  0.  (16) 

^From  the  partial  differential  equation  point  of  view,  this  might  well  seem  the  correct  result; 
V  •  B  =  0  is  an  initial  condition,  and  this  equation  guarantees  that  V  •  B  =  0  is  maintained 
throughout  the  evolution.  For  the  modified  form  of  the  equations,  the  evolution  equation 
for  the  magnetic  field  is 


VB  -b  BV  u  -  B  Vu 


)  = 


—  (V-B)-bV-(uV-B)  =  0. 

at 
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Thus  the  addition  of  the  source  term  has  modified  the  evolution  equation  for  V  •  B  so  that 
the  quantity  V  •  B//»  is  treated  as  a  passive  scalar.  This  is  clearly  the  more  numerically 
stable  of  the  two  evolution  equations;  any  local  V  ■  B  that  is  created  is  convected  away. 

The  above  derivation  gives  all  the  pieces  for  building  an  ideal  MHD  solver  that  works 
for  two-dimensional  problems,  without  having  to  resort  to  non-collocated  variables  or  a 
projection  algorithm.  Specifically,  a  Floe-type  approximate  Riemann  solver  has  been  imple¬ 
mented,  where  the  wave  strengths  and  speeds  are  derived  from  the  above  left  eigenvectors 
and  eigenvalues.  The  eigenvectors  are  properly  normalized  to  avoid  difficulties  associated 
with  coinciding  wave  speeds  [5].  The  average  state  needed  at  cell  interfaces  is  computed  by  a 
simple  average  of  left  and  right  states  (although  a  Roe  average  does  exist  for  the  ideal  MHD 
equations  [4]).  The  source  term,  though  small,  is  calculated  in  each  cell,  and  added  to  the 
residual.  The  resulting  code  is  first  order  in  space  and  time. 

3  A  Test  of  the  Eight- Wave  Riemann  Solver 

Brio  and  Wu  [2]  developed  a  test  problem  for  one-dimensional  MHD  solvers  based  on  the 
shock-tube  problem  of  Sod  [6].  Two  stationary  plasmas  are  separated  by  a  membrane  which 
is  removed  at  <  =  0,  allowing  the  plasmas  to  interact.  The  test  problem  used  here  for  the 
two-dimensional  MHD  solver  is  a  rotated  version  of  the  Brio-Wu  problem.  The  left  and 
right  input  states,  and  the  orientation  of  propagation  of  disturbances  to  the  grid,  is  shown 
in  Figure  2.  In  the  Brio-Wu  problem  (the  top  figure),  the  boundary  conditions  are  that  the 
problem  is  periodic  along  a  line  y  =  constant;  in  the  current  test  problem  (the  bottom  figure), 
the  boundary  conditions  are  that  the  problem  is  periodic  along  a  line  x  -{-y  =  constant. 

Both  the  rotated  and  non-rotated  problems  were  run  on  coarse  (600  cells  in  x)  and  fine 
(1200  cells  in  x)  grids.  The  time  step  was  taken  as  AtfAx  =  0.2,  which  corresponds  to  a 
CFL  number  of  approximately  0.8  on  the  non-rotated  problem.  The  ratio  of  specific  heats, 
7,  was  2.0.  The  number  of  time  steps  taken  on  the  coarse  and  fine  grids  were  100  and  200, 
respectively.  The  a:— axis  in  the  plotted  results  from  the  rotated  problem  was  scaled  by  a 
factor  of  V2,  to  account  for  the  fact  that  the  CFL  number  is  lower  for  the  rotated  problem 
than  for  the  non-rotated  problem. 

Figures  .3-7  show  comparisons  of  the  results  on  the  fine  grid  of  the  non-rotated  shock-tube 
problem  with  the  (scaled)  results  of  the  rotated  shock-tube  problem  for 

3.  density  (/?); 

4.  pressure  (p); 

5.  velocity  component  normal  to  the  original  discontinuity  (««); 
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p=l  pssl  p=l/8  p=l/10 


B  =3/4  B  =1  B  =3/4  B  =-l 
X  y  X  y 


p=l  p=l  p=l/8  p=l/10 


B— *1  g-.7  g_7  g_*l 

*  4V2  y  4V2  *  4V2  y  4V2 


Figure  2:  A  Test  Problem  for  Two-Dimensional  MHD 
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Figure  3:  Density  in  the  Rotated  and  Non-Rotated  Shock  Tubes 
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Figure  4:  Pressure  in  the  Rotated  and  Non-Rotated  Shock  Tubes 
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Figure  5:  Normal  Velocity  in  the  Rotated  and  Non-Rotated  Shock  Tubes 
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_ non-rototed 

_ rotated 
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Figure  6:  Tangential  Velocity  in  the  Rotated  and  Non-Rotated  Shock  Tubes 


X 


Figure  7:  Tangential  Magnetic  Field  in  the  Rotated  and  Non-Rotated  Shock  Tubes 
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Figure  8:  Normal  Magnetic  Field  in  the  Rotated  Shock  Tube  (Coarse  and  Fine) 


6.  velocity  component  tangential  to  the  original  discontinuity  (ut); 

7.  magnetic-held  component  tangential  to  the  original  discontinuity  {Bt). 

As  can  be  seen,  the  agreement  is  quite  good,  with  the  results  of  the  two  cases  nearly  indis¬ 
tinguishable  for  all  but  the  normal  component  of  velocity.  The  errors  in  u„  are  balanced  by 
errors  in  the  magnetic-field  component  normal  to  the  original  discontinuity  (B„).  Figure  8 
shows  Bn  for  the  rotated  shock-tube  problem  on  the  coarse  and  fine  grids.  In  the  non-rotated 
problem,  B„  =  0.75  throughout  the  tube.  As  can  be  seen,  there  are  errors  on  the  order  of  a 
few  percent  in  Bn  on  the  coarse  grid,  but  the  errors  are  reduced  as  the  grid  is  refined. 

4  Concluding  Remarks 

In  some  respects,  this  paper  presents  the  development  of  only  one-eighth  of  a  Riemann  solver. 
Seven  of  the  eight  waves  of  the  Riemann  solver  are  the  same  as  those  used  in  previous  work  on 
upwind  methods  for  MHD.  The  deceptively  simple  eighth  wave  that  arises  from  the  analysis, 
however,  is  of  a  different  chfiracter  than  the  other  seven  —  it  arises  only  in  multi-dimensional 
problems,  and  it  is  crucial  for  understanding  and  solving  those  problems.  It  plays  the  very 
important  role  of  stabilizing  the  numerical  method  with  respect  to  the  small  amounts  of 
V  ■  B  generated  in  solving  the  discrete  MHD  equations. 

Given  the  meteoric  rise  of  Riemann  solvers  in  the  computation  of  compressible  gas  dy¬ 
namics,  it  is  not  very  risky  to  predict  that  schemes  based  on  Riemann  solvers  will  play  an 
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increasingly  important  role  in  the  computation  of  compressible  conducting  flows.  The  ability 
of  Riemann  solvers  to  capture  strong  discontinuities  robustly  and  with  minimal  dissipation, 
the  framework  that  Riemann  solvers  provide  for  implementing  stable  boundary  procedures, 
and  the  aesthetically  appealing  physical  basis  of  Riemann  solvers  are  all  strong  arguments 
for  their  use.  The  aim  of  this  paper  is  to  remove  what  is  hopefully  one  of  the  last  major 
obstacles  to  the  use  of  Riemann  solvers  in  large-scale  codes  for  computing  multi-dimensional 
conducting  flows. 
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